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ABSTRACT 

We show that the gas giant exoplanet HD 189733b is less oblate than Saturn, based on Spitzer Space Tele- 
scope photometry of seven transits. The observable manifestations of oblateness would have been slight anoma- 
lies during the ingress and egress phases, as well as variations in the transit depth due to spin precession. Our 
nondetection of these effects gives the first empirical constraints on the shape of an exoplanet. The results are 
consistent with the theoretical expectation that the planetary rotation period and orbital period are synchro- 
nized, in which case the oblateness would be an order of magnitude smaller than our upper limits. Conversely, 
if HD 189733b is assumed to be in a synchronous, zero-obliquity state, then the data give an upper bound 
on the quadrupole moment of the planet (J2 < 0.068 with 95% confidence) that is too weak to constrain the 
interior structure of the planet. An Appendix describes a fast algorithm for computing the transit light curve of 
an oblate planet, which was necessary for our analysis. 

Subject headings: methods: numerical — stars: planetary systems — techniques: photometric 



1. INTRODUCTION 

Planets are not exactly spherical. Departures from spheric- 
ity are caused by rotation, external gravitational tides, and 
material rigidity. For gas giant planets, material rigidity is 
not important, but tides and rotation are very important. The 
tidal bulges raised on "hot Jupiters" by their parent stars al- 
low angular momentum to be exchanged between the plan- 
etary rotational and orbital motions, and result in long-term 
energy dissipation (Goldreich & Soter 1966, Peale 1999, Mur- 
ray & Dermott 2000). As for rotation, both Jupiter and Saturn 
are visibly flattened due to centrifugal forces, with their polar 
radii shorter than equatorial radii by 6.5% and 9.8%, respec- 
tively (Lindal et al. 1981, 1985). 

It would be of great interest to determine the shapes of plan- 
ets outside of the Solar System. Knowledge of the equilibrium 
shape of an exoplanet would provide information on its rota- 
tion rate and internal density structure, which in turn would 
give clues about its formation and evolution. One way this 
knowledge might be gained is through precise photometry of 
exoplanetary transits. Seager & Hui (2002) and Barnes & 
Fortney (2003) calculated the difference between the transit 
light curve of an oblate planet and of a spherical planet with 
the same cross-sectional area. Unfortunately the oblateness- 
induced signal is expected to be quite small for most of the 
currently known transiting planets because their rotation peri- 
ods are expected to be tidally synchronized with their ^3 day 
orbital periods (much longer than the «10 hr rotation peri- 
ods of Jupiter and Saturn). For the representative case of 
HD 209458b, Barnes & Fortney (2003) found the amplitude 
of the theoretically expected oblateness-induced signal to be 
0.1 /xmag, well below the limiting precision of any current 
photometer. 

The goal of this work was to place the first empirical con- 
straints on the shape of an exoplanet. While the assumption 
of spin-orbit synchronization seems reasonable, it is worth 
checking on such assumptions whenever possible. It is also 
useful to know the upper limits that can be achieved with 
current data, with an eye toward planning future observations 
and photometric instruments. By selecting the most favorable 



planet and the best available data, we find that an oblateness 
as large as that of Saturn can be ruled out, and an oblateness 
as large as that of Jupiter is also ruled out if it is accompanied 
by a moderate obliquity. In addition, we model the effects of 
spin precession on transit light curves, which were not con- 
sidered by Seager & Hui (2002) or Barnes & Fortney (2003), 
and provide a fast algorithm for calculating the light curve of 
an oblate planet. 

This paper is organized as follows. In §[2] we review the 
relevant physics and geometry. In § [3] we present seven 
Spitzer 8 fim transit observations of the well-studied system 
HD 189733 (Bouchy et al. 2005). In §H we present the re- 
sults of fitting a parameterized model to the data, including the 
sky -projected oblateness and obliquity in addition to the usual 
transit parameters. We consider the case of a fixed orientation 
for the plan et (§ 14. 11 ) as well as uniform precession of the spin 
axis (§ 14.21 ). In the latter case we are able to constrain the true 
oblateness and obliquity as well as the precession period. In 
§[5] we summarize our methods and results and suggest some 
possible extensions. Appendix lAl describes the algorithm we 
used to produce transit light curves of oblate planets. 

2. REVIEW OF PLANETARY OBLATENESS 
2.1. Rotation 

A uniformly rotating self-gravitating fluid takes on the fig- 
ure of an oblate spheroid, with its minimum diameter along 
the axis of rotation (Eddington 1926). The shape may be 
quantified by the oblateness (or flattening) parameter /, de- 
fined as 
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where R eq and R po \ are the equatorial and polar radii, respec- 
tively (Murray & Dermott 2000). The angle between the polar 
axis and the orbital axis is the obliquity 6. 

The relationship between / and the rotation period P rot is 
found by identifying the contours of constant potential, where 
the potential has both gravitational and centrifugal terms. The 
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gravitational terms are 
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where J„ are the spherical mass moments 1 associated with ro- 
tation, if) is the planetary colatitude, and V„ is the Legendre 
polynomial of degree n. The centrifugal terms are 



V c (r,V)= -tt 2 r 2 [7> 2 (cosV)-l], 



(3) 



where fl = 2n/P ml is the rotational angular frequency. As- 
suming that the quadrupole moment Ji is the most important 
moment, the total potential is 
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Because the surface of the planet lies on an equipotential, 
V to t(fle q ,§) = Vto t (flpoi,0), giving 
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This equality implies (to leading order in f) 
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Derivations of this relation are also given by Murray & Der- 
mott (2000) and Hubbard (1984). Solving for P rot = 2tt/U we 
find 
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This result reduces to Eqn. (6) of Seager & Hui (2002) if the 
J2 term is neglected. 

For Solar system planets, / is measured from direct images 
and J2 is measured by monitoring elliptical orbits of satellites 
whose orbits precess in response to the aspherical gravita- 
tional field. Table[T]gives / and J 2 for planets in our Solar Sys- 
tem, to help place our results for the exoplanet HD 1 89733b in 
context. The maximum possible value for / occurs at the ro- 
tational breakup limit, when the outward centrifugal accelera- 
tion equals the gravitational acceleration at the equator. With 
reference to Eqns. (fJlO, this criterion gives / < 5 + 1/2- 

Another relation between J 2 and / can be obtained by as- 
suming the planet's interior to be in a state of hydrostatic equi- 
librium, leading to the Darwin-Radau approximation, 
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Here, C is the dimensionless moment of inertia, defined such 
that <CM p rL is the moment of inertia about the spin axis 
(Murray & Dermott 2000). The moment of inertia of Solar 

' J " = So* /-l r"V„(v)p(r,n)27rr 2 dy. dr 



TABLE 1 
Shape parameters of Solar 
system planets 



Planet 


Oblateness / 


h 


Mercury 


0.00012 


0.000060 


Venus 


0.00009 


0.000004 


Earth 


0.00350 


0.001083 


Mars 


0.00520 


0.001960 


Jupiter 


0.06487 


0.014736 


Saturn 


0.09796 


0.016298 


Uranus 


0.02293 


0.003343 


Neptune 


0.01708 


0.00341 1 



REFERENCES. — Murray & Der- 
mott (2000), Barnes & Fortney (2003), 
Hubbard (1984). 

system planets is not directly measurable but models of gas 
giant interiors suggest C sw 0.23 (Hubbard & Marley 1989). 

2.2. Sky projection 

A transit by an oblate planet will produce a slightly differ- 
ent light curve than a transit by a spherical planet of the same 
cross-sectional area. The differences are most pronounced 
during the ingress and egress phases of the transit. There are 
also slight differences during the complete phase of the tran- 
sit, due to the limb darkening of the stellar photosphere. For 
illustrations of these effects we refer the reader to Seager & 
Hui (2002) and Barnes & Fortney (2003). 

The shape of a single transit light curve will depend only 
on the sky projection of the planet's figure at the time of the 
transit. The sky projection of an oblate or prolate spheroid (or 
any triaxial ellipsoid) is an ellipse. We define the projected 
oblateness fx as (a - b) /a, where a and b are the lengths of 
the major and minor axes of the ellipse. We also define the 
projected obliquity 9± as the angle between the major axis 
and the transit chord. The relations between these projected 
quantities and the unprojected quantities / and 9 can be found 
with elementary geometry: 



/ ± = l-ysin 2 c?' + (l-/) 2 cos 2 (?', (9) 

tan6>j_ =tanf? sin</>, (10) 

where 9' is the angle between the planetary rotation axis and 
the sky plane, given by 



cos 2 9' = sin 2 9 sin 2 6+ cos 2 ( 
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and (j> is the azimuthal angle of the line of nodes between 
the planetary equator and the orbital plane, with <fi = cor- 
responding to the case when the rotation axis is tipped toward 
the observer. A derivation of Eqn. (0 is also given by Barnes 
(2009). 

The calculation of the transit light curve of an oblate planet 
is computationally intensive, because the intersection points 
between a circle and an ellipse are nonanalytic, and because 
the precision of the calculation must be very high in order to 
isolate the small oblateness-specific effects. In Appendix lAl 
we describe a fast algorithm that we developed for our study, 
building on the previous work by Seager & Hui (2002) and 
Barnes & Fortney (2003). 

2.3. Spin precession 

External gravitational forces will cause a planet's spin axis 
to precess. In particular, the spin axis of a rotationally- 
induced, oblate spheroidal planet will precess in response to 
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the gravitational torque of its host star, with a period 

p = 2P °rbC 1 

prec 3 P rot J 2 cos d 



(12) 



where 9 is the planetary obliquity (Ward 1975). This expres- 
sion assumes that the orbit is fixed, and in particular that the 
orbital angular momentum is much larger than the rotational 
angular momentum, a good approximation even for close-in 
giant planets rotating as fast as Jupiter or Saturn. 2 

For Solar system planets, the precession periods are much 
longer than 1 yr. For example, Saturn completes one preces- 
sion cycle in w 7 x 10 6 yr (ignoring the effect of moons and 
rings; Ward & Hamilon 2004). However, because P prec scales 
as P^ rb , sufficiently close-in planets would precess with pe- 
riods short enough to be directly observable. If Saturn's or- 
bital period were chan ged t o P or b = 3 days it would precess 
with P prec ~ 0.6 yr. In § 14.21 we consider how spin precession 
would be manifested in a collection of transit light curves. 

2.4. Tidal deformation 

Hot Jupiters are stretched radially (along the radius vector 
of the orbit) due to the gravitational tide from the nearby host 
star. Tidal dissipation causes important long-term changes in 
the orbital parameters, including spin-orbit synchronization 
and orbital circularization. The characteristic timescale for 
synchronization is 
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where luq is the initial angular rotation frequency, M p is the 
planet's mass, M t is the stellar mass, and Q is the specific dis- 
sipation factor of the tidal oscillator (Goldreich & Soter 1966, 
Hubbard 1984, Guillot et al. 1996, Murray & Dermott 2000). 
Many investigators adopt Q ~ 10 5-6 for gas giant planets, al- 
though there is little empirical information about Q and some 
indications that Q may be larger (see, e.g., Hellier et al. 2009). 

The obliquity is also expected to be driven to zero as a result 
of tidal evolution, although in an interesting and nonuniform 
manner (Peale 1999). It is possible to maintain a nonzero 
obliquity in a so-called Cassini state (Ward 1975, Winn & 
Holman 2005) but for hot Jupiters the high-obliquity states 
are unstable (Fabrycky et al. 2007, Levrard et al. 2007). 

It would be difficult to detect a radially-directed tidal bulge 
through transit photometry. For tidally synchronized planets, 
the amplitude of the rotational deformation is roughly 1/3 
that of the tidal component (Murray & Dermott 2000, pg. 
156). However, the long axis is closely aligned with the line 
of sight during transits, and consequently the sky projection 
of the tidal bulge is smaller than the bulge itself by a factor 
a/R+. This causes a reduction in the transit signal of the ra- 
dial bulge by an order of magnitude or more, compared to 
the rotational bulge. Therefore in this paper we consider only 
the deformation associated with rotation. Ragozzine & Wolf 
(2009) describe some potentially observable consequences of 
the radial bulge. 

2 In reality the orbital and spin axes both precess about the total angular 
momentum vector. The nodal precession of the orbit would be detectable in 
principle through changes in the transit impact parameter. However, even if 
HD 189733b were spinning as fast as Jupiter, the orbital angular momentum 
would exceed the rotational angular momentum by a factor of 1000, and the 
resulting inclination variation would not be detectable in the current data. 



2.5. Expectations for HD 189733b 

In this section we use the preceding formalism and the algo- 
rithm described in the Appendix to compute theoretical light 
curves for the transiting exoplanet HD 189733b. This par- 
ticular exoplanet was chosen because it is the most favorable 
case currently known for the detection of oblateness-induced 
signatures, owing to the bright parent star, large transit depth, 
midrange transit impact parameter, and the availability of a 
large corpus of high-precision transit data. We took the sys- 
tem parameters from Torres et al. (2008) and considered two 
different cases for the oblateness and obliquity. 

First, we imagined that HD 189733b is as oblate as Saturn 
(f± = 0.098), with a projected obliquity of 9± = 45°. To iso- 
late the signal due to oblateness, we calculated the theoretical 
transit light curve for the oblate planet, and then subtracted 
the best-fitting model for a spherical planet, following Barnes 
& Fortney (2003). The residuals are shown in the left panel 
of Figure Q] An error bar is also shown, representing the fore- 
casted precision of a 2 min sample based on the Spitzer data 
that are currently available for this planet (see § 3). This cal- 
culation suggested that the Spitzer data would be capable of 
placing physically meaningful constraints on the oblateness, 
and motivated our further study. 

Second, we assumed that HD 189733b has been tidally 
synchronized with P rot = f or b = 2.2 days. This is reason- 
able because for Q = 10 5 and uj pI i m = 1.7 x 10~ 4 s" 1 (Jupiter's 
estimated values; Guillot et al. 1996) we find r ~ 10 6 yr, 
much shorter than the estimated few-Gyr main-sequence age 
of the star (Torres et al. 2008). Furthermore, the observa- 
tion that the orbit of HD 189733b is nearly circular (Winn et 
al. 2007b, Knutson et al. 2007a) is independent evidence for 
tidal evolution, and the theoretical synchronization timescale 
is shorter than the circularization timescale. We estimated the 
rotationally-induced oblateness using the Darwin-Radau re- 
lation and Eqn. ©, finding / ps 0.003. The right panel of 
Figure Q] shows the difference between a light curve of an 
oblate planet with f± = 0.003 and 9± = 0, and the best-fitting 
model for a spherical planet. The peak-to-peak amplitude of 
the residuals is approximately 2 x 10~ 6 , below the precision 
of the available data. Thus, if spin-orbit synchronization has 
been achieved, we would expect a null result from an analysis 
of the Spitzer data. 

3. OBSERVATIONS AND DATA REDUCTION 

Transits of HD 189733b have been observed with many 
instruments and in many bandpasses (see, e.g., Bouchy et 
al. 2005, Winn et al. 2007b, Pont et al. 2007, Beaulieu et 
al. 2008, Miller-Ricci et al. 2008, Desert et al. 2009). For our 
study the best available data are the 7 different transits that 
were observed with the Spitzer Space Telescope at a wave- 
length of 8 /im, with the InfraRed Array Camera (IRAC; 
Fazio et al 2004). They form a homogeneous and precise data 
set, and by virtue of the relatively long observing wavelength 
the data are less affected by stellar variability, star spots and 
limb darkening, which would otherwise confound the attempt 
to detect or constrain oblateness. The first of the 7 IRAC time 
series was presented by Knutson et al. (2007a), and the other 
6 time series are based on observations by Agol et al. (2009). 
The elapsed time between the first and seventh transits was 
1.6 yr or 268 transits. Assigning epoch zero to the second 
transit (on UT 30 June 2007) the observed transits have epoch 
numbers -110,0, 1,51,62, 157, and 158. 

All 7 transits were observed with the 8 fim channel of the 
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FIG. 1. — The theoretical oblateness-induced signal in the transit light curve of HD 189733b. Each panel shows the difference in the time series between the 
theoretical transit light curve of an oblate planet, and the best-fitting light curve of a spherical planet. In computing the theoretical light curves, all the system 
parameters except oblateness and obliquity were taken from Torres et al. (2008). Left. — Assuming HD 189733b is as oblate as Saturn and is more oblique 
(fj_ = 0.098, 6± = 45°). The error bar indicates the expected Icr noise level of a composite light curve built from seven Spitzer observations, binned into 2 min 
intervals. Right. — Assuming f± = 0.003 and 6± = 0, as appropriate for the case of a tidally synchronized planet. 



IRAC instrument in subarray mode, in which only a 32x32 
pixel subarray of the detector is recorded. Images are obtained 
with a cadence of 0.40 s (with an integration time of 0.32 s). 
The data are packaged into post-calibration FITS files, each 
consisting of 64 images (representing a total of 26.5 s of in- 
tegration time). Approximately 500 files span the transit of 
HD 189733b. In addition to providing excellent time reso- 
lution, the subarray mode avoids saturation even for stars as 
bright as HD 189733 (V = 7.7, H0g et al. 2000). 

We began our data reduction with the post-calibration im- 
ages, downloaded from the Spitzer data archive. First, we in- 
spected each block of 64 images and determined the centroid 
of the target star in each image. We excluded any image in 
which the centroid differed by more than 0.05 pixels from the 
median position of all 64 images. Fewer than 1% of the im- 
ages were excluded by this criterion, except for epochs 1 and 
157 in which 8% and 5% of images were excluded, respec- 
tively. Next, we formed a mean image based on each block of 
64 images, disregarding any pixels whose values deviated by 
more than 3.5cr from the median value for that pixel. Fewer 
than 1 % of the pixels were excluded by this criterion. 

We then performed aperture photometry on the mean im- 
ages, using a circular aperture of radius 4.5 pixels centered 
on HD 189733. To estimate the background level, we also 
summed the flux within several rectangular apertures located 
far away from both HD 189733 and its fainter companion. 
The background level was subtracted from the aperture sum. 
At this point, we had a time sequence of measurements of the 
relative flux density of HD 189733. 

Among users of the IRAC 8 /im channel there is a well- 
known systematic effect which has been called the "ramp," 
because it is manifested as a gradual rise in the count rate at 
the beginning of an observation. It is attributed to charge trap- 
ping in the detector. It is generally modeled as a multiplica- 
tive, time-variable correction (see, e.g., Knutson et al. 2007b, 
Gillon et al. 2007, Nutzman et al. 2008). For each data set, 
we clipped the most strongly-varying portion of the ram p, and 
modeled the rest as a quadratic function of time (see § 14.1b . 



We also clipped the long post-egress portion of the time se- 
ries by Knutson et al. (2007a). Figure [2] shows the final time 
series, after correcting for the ramp. 

We assessed the noise characteristics of each time series 
with the solveredwv IDL routine 3 . This algorithm fits 
time-series data with a model in which the noise is an ad- 
ditive combination of white noise and 1 // noise. The ampli- 
tudes of each noise component are estimated from the data as 
described by Carter & Winn (2009). Many of our reduction 
parameters (aperture size, thresholds for clipping, etc.) were 
chosen by attempting to minimize the time-correlated com- 
ponent of the noise. Our final time series had white-noise 
amplitudes of 586, 571, 642, 544, 560, 607 and 536 ppm, 
for epochs -110, 0, 1, 51, 62, 157, and 158, respectively. 
In all cases the amplitude of the 1// component was smaller 
than 2 ppm. The theoretical limiting precision due to photon- 
counting noise was approximately 460 ppm. 

4. ANALYSIS 
4. 1 . Fixed orientation 

In the analysis described in this section, we assumed the 
orientation of HD 189733b to be fixed in space over the 1.6 yr 
span of the observations. The projected oblateness parameters 
f± and 6± were taken to be constants. In the next section we 
will describe an analysis in which that restriction was lifted. 

The first step was a preliminary fit to all of the data, with the 
goal of determining the midtransit time, out-of-transit flux, 
and the coefficients of the ramp correction function for each 
epoch. These quantities are weakly correlated with the other 
parameters describing the light curves (including the pro- 
jected oblateness and obliquity), and thus may be fixed in the 
subsequent analysis without significantly affecting the results. 

The ramp correction function was 

C ramp (f;c ,ci;fo)= 1 +c (f-f ) + ci(f-f ) 2 (14) 
where cq and c\ are adjustable parameters, and frj is a partic- 

3 |http: //www.mit ■ edu/~earter ja/cocle/| 
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ular time near midtransit. In addition to the ramp correction, 
each light curve was described by an out-of-transit flux level 
Fq and midtransit time fy. There were also 7 parameters com- 
mon to all the light curves: the mean projected radius ratio 
R±/R+ (where R± = R^^l —f±), orbital inclination i, nor- 
malized orbital distance a/R±, quadratic limb-darkening co- 
efficients Mi and U2, projected oblateness f±, and projected 
obliquity 9±. The limb darkening coefficients parameterize 
the stellar brightness profile, /*(r; Wi, U2), defined as 

4(/i;Mi,M 2 ) = 4(l) [l-M 1 (l-/i)-M 2 (l-M) 2 ] (15) 



where fi = Vl — r 2 and r is the sky-projected distance from the 
center of the star. 

To derive the best-fitting parameter values we minimized 
the standard \ 2 statistic, 
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where o indexes the observation number, N„ is the number of 
data points in observation o, f)°(obs) is a measurement of the 
relative flux of HD 189733 during observation o, fy'(calc) is 
the calculated flux at that time according to the model, and 
a" is the uncertainty in flux measurement during observation 
o. We took a° to be the white noise amplitude of each time 
series as specified in § [3] As described above, we found the 
correlated noise component to be negligible. We minimized 
X 2 using the AMOEBA routine (Press et al. 2007). The ramp- 
corrected light curves and the best-fitting model light curves 
are shown in Figure [2] 

Then, after fixing the midtransit times, ramp correction pa- 
rameters, and out-of-transit flux levels at the best-fitting val- 
ues, we determined the credible intervals for the other pa- 
rameters using a Markov Chain Monte Carlo (MCMC) algo- 
rithm. 4 We found that the quadratic limb darkening parameter 
was very poorly constrained by the data, and therefore for sub- 
sequent work we fixed that parameter at the value obtained in 
the preliminary fit. The parameters that were allowed to vary 
in the chain were R±/R+, i, a/R+, u\, Tq, f±, and 9±. Here, 
7b represents an overall timing offset of all 7 transits, which 
is needed because the oblateness parameter is correlated with 
an overall time shift. 

In the MCMC, the jump-transition probability was propor- 
tional to exp(~x 2 /2) with % 2 given in Eqn. (fTol l. We used 
Gibbs sampling in the construction of a chain of 5 x 10 5 links. 
We selected individual parameter jump sizes such that the 
fraction of jumps accepted by a Metropolis-Hasting condition 
was between 30% and 60% for each parameter. Uniform pri- 
ors were used for all parameters. We verified that the resulting 
posterior distributions had converged sufficiently. 

Table[2]gives the results. For most parameters we report the 
median value of the posterior probability distribution, along 
with error bars defined by the 15.85% and 84.15% levels of 
the cumulative distribution. For the projected oblateness, we 
report the 95% -confidence upper limit. The projected obliq- 
uity was unconstrained. Figure (f3]l shows some of the poste- 
rior probability distributions. Figure |4] shows the constraints 
in the f±-9± plane. Only the results for 9± > are shown, as 
the constraints were found to be symmetric about 9± = 0. For 

4 For background on the MCMC method, see Gregory (2005), and for 
examples of applications to transit light curves, see Holman et al. (2006), 
Winn et al. (2007a), or Burke et al. (2007). 



comparison, we have indicated in the figure the oblatenesses 
of Saturn, Jupiter and Uranus. 

As anticipated from the theoretical work of Seager & 
Hui (2002) and Barnes & Fortney (2003), the projected 
oblateness is most tightly constrained when the projected 
obliquity is near 45°. We determined the 68% and 95% con- 
fidence regions in the f±-9± plane with the following tech- 
nique. First, we divided the projected obliquity range into 
bins of constant width. Next, for the chain links in each bin, 
we sorted the values of f± and determined f' ± such that 68% 
(95%) of the values were less than Fig. (0]i shows the 
resulting confidence curves. 

Our analysis rules out a projected oblateness that is equal to 
Saturn's oblateness with >95% confidence, regardless of the 
projected obliquity. A projected oblateness equal to Jupiter's 
oblateness is also ruled out with >95% confidence, except for 
obliquities within about 7° of either 0° or 90°. For a pro- 
jected obliquity near 45°, we may also exclude a projected 
oblateness equal to the oblateness of Uranus, with 95% con- 
fidence. Our constraints are consistent with the theoretical 
expectation that HD 189733b is rotating in synchrony with its 
orbit (f ± fa / ps 3 x 10" 3 and 9j_ m 0°). 

The nondetection of oblateness suggests that the planet can- 
not be rotating too quickly, but the results cannot be translated 
directly into a lower bound on P TOt because of the sky projec- 
tion. Even a very rapidly rotating (and very oblate) planet is 
consistent with the data as long as the planet's rotation axis 
is pointing at the observer (9 = 90° and = 0), leading to a 
circular projected figure. If we assume 9 ps then we may set 
an upper bound on P mt . We further assume that J2 is given by 
the Darwin-Radau relation (Eqn. [HJ, with C = 0.225. Under 
these assumptions P mt > 0.39 d (9.4 hr) with 95% confidence. 
Larger values of C would correspond to longer rotational pe- 
riods. The median value of the posterior distribution is P mt 
is 0.65 days, but we do not attribute any significance to that 
result, as the shape of the posterior distribution (including the 
median) is strongly affected by our assumption of a uniform 
prior in f± . The lower bound on P m is less sensitive to the 
prior. 

Conversely, if we assume HD 189733b to be synchronously 
rotating with zero obliquity and P rot = 2.218573 days, then we 
may place an upper bound on the rotationally-induced J2, us- 
ing Eqn. ©. The resulting posterior probability distribution 
for J2 is shown in Figure We find that Jo must be smaller 
than 0.068 with 95% confidence. To place this in perspective, 
we note that if HD 189733b were a uniform-density sphere 
then under the same assumption of spin-synchrony and zero 
obliquity one would expect Jo = 0.0018. 5 More centrally con- 
densed planets would lead to smaller values of J%. Hence, the 
empirical upper bound on J2 is not constraining on physically 
plausible models of the planet's interior. 

4.2. Uniform precession 

As discussed in § 12.31 if HD 1 89733b is significantly oblate 
and oblique, then the spin axis will precess with a much 
shorter period than the precession periods of Saturn or Jupiter. 
The previous analysis assumed a fixed orientation for the 
planetary spin axis, and therefore the results must be under- 
stood as valid only for cases in which the precession period is 
much longer than 1 yr. The actual precession period depends 
on the internal constitution of the planet, specifically its mo- 

5 This was derived from Eqn. (3) of Ragozzine & Wolf (2009), using a 
Love number k% = 1.5 for a uniform-density sphere. 
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FIG. 2. — The relative 8 /im brightness of HD 189733 during 7 different transits of its giant planet, as measured w ith S pitzerflRAC. The time series were 
corrected for the detector "ramp," and offset in flux for display purposes. The solid curve is the best-fitting model. See jj |4. ll for details. 



ment of inertia and J 2 ■ In the next part of our analysis we 
relaxed the assumption of a fixed orientation for the planetary 
spin axis, and allowed the axis to precess uniformly about the 
orbital axis. 

We assumed that the true oblateness / and true obliquity 9 
are constant in time, and allowed the azimuthal angle <fi to be 
a linear function of time, 

2lTt 

m=- r - +00, (17) 

*prec 

for some precession period P prec and initial phase angle 4>q. 
The time evolution of fx and 6± are then given by Eqns. (|9} 
ITOb . A consequence of this time evolution is that the transit 
depth, 6(t), is variable. Apart from small corrections due to 
limb darkening, the transit depth is equal to the ratio of areas 
of the projected disks of the planet and star: 

m ~(lt) [ 1_/±(f) ]' (18) 



For large oblateness and obliquity, the transit depth variations 
can be easier to detect than the slight distortions in the ingress 
and egress portions of a single light curve. For a Saturn-like 
oblateness (/ as 0.1) at 45° obliquity, the transit depth vari- 
ations are nearly 5% (see Fig.|6]l. For the 7 Spitzer time se- 
ries, the transit depths were found to agree with one another 
to within 0.5%. This suggests that if HD 189733b were as 
oblate as Saturn and were also significantly oblique, then the 
spin-precession period must be much longer than 1 .6 yr. For a 
quantitative and physically self-consistent analysis, we com- 
puted a second MCMC, using the parameters /, 8, P prec and 
4>q to describe the shape of the planet, 6 rather than the sky 
projected parameters f± and 9±. 

In principle, we could have used Eqn. ( fl2] i to enforce a re- 
lationship between P prec and /, but that would have required 
some assumptions about the interior density structure of the 
planet. We preferred to allow P prec to be a free parameter 

6 It proved advantageous to use i? c q V 1 - f as a fitting parameter, rather 
than R c q, to reduce the correlation with /. 
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TABLE 2 

Model Results: Static Shape Parameters 



Parameter Median Uncertainty 



Transit parameters: 

Rj_/R i , 0.154679 

Orbital Inclination, i [degrees] 85.749 

a/R* " 8.924 

Projected oblateness /j_ 

Projected obliquity 9± unconstrained 



±0.000067 

±0.026 

±0.022 

< 0.056 (95% conf.) 



Limb darkening parameter «i 
Limb darkening parameter uj 
Midtransit time shift Tq [s] 



0.076 ±0.011 

0.034 (fixed) 

0.0 ±1.3 



Derived parameters" : 
Rotational Period [days]*- c - > 0.39 (95% conf.) 

h A - < 6.8 x 10~ 2 (95% conf.) 



NOTE. — (a) Assuming zero obliquity, (b) Calculated using Eqn. {7) with 
i?* = 0.756 if 0, Mp = 1.144 M Jup (Torres et al. 2008). (c) Assuming the validity 
of the Darwin-Radau approximation (Eqn. [8) with C = 0.225. (d) Calculated 
using Eqn. (5) assuming P rot = P OI t = 2.218573 days. 





i [degrees] Projected-obliquity 8 1 [degrees] 



FIG. 3. — Joint confidence regions and posterior distributions for transit parameters and the projected shape parameters of HD 189733b, based on Spitzer 
observations of 7 transits. The shape parameters were assumed to be time-independent. 
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FIG. 4. — Posterior probability distribution for the projected oblateness and obliquity of HD 189733b, based on Spitzer observations of 7 transits. The shape 
parameters were assumed to be time-independent. The solid blue curves bound the regions containing 68% or 95% of the probability in the (f±—9±) plane, 
marginalized over all other parameters. The black points are 10,000 representative values from the Markov chain, to illustrate the probability density. The star 
shows the expected shape parameters if HD 189733b is spin-orbit synchronized. The dashed lines indicate the oblateness of Jupiter, Saturn, and Uranus, for 
comparison. 
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FIG. 5. — Posterior probability distribution for the quadrupole moment (J2) 
of HD 189733b. See §|4]for a discussion of the underlying assumptions in 
this calculation. 

to avoid making such assumptions. The structure of close- 
in planets may differ from that of Jupiter and Saturn, due to 
a different history of formation and evolution, as well as the 
tidal influence of the parent star. We restricted our attention to 
^prec < 10 yr, because for longer precession periods the orien- 
tation of the spin axis would appear fixed over the 1.6 yr span 
of the Spitzer observations, and therefore our analysis would 
revert to that of the previous section. 

We created two Markov chains, each with 10 6 links, and 
merged them after removing the first 10 5 links of each chain. 
The results are depicted in Fig.[6]as confidence regions in the 
f-9 plane. This figure is analogous to Fig. @, although in 
this case we are constraining the true oblateness and obliquity 
as opposed to the projected quantities. Solutions with high 
oblateness and obliquity are ruled out by the observed ab- 
sence of transit depth variations. Indeed, for f prec < 10 yr we 
can rule out a wider swathe of obliquity-oblateness parameter 
space than we did under the assumption of a nonprecessing 
planet. 

For reference, Fig. [6] also shows contours indicating the 
amplitude of the expected transit depth variations, as well 
as contours indicating the spin precession period calculated 
from Eqn. il2i , assuming C = 0.225. The precession period 
is smaller than 10 yr over more than 90% of the f-9 plane. 

5. DISCUSSION 

We have made the first attempt to measure the shape of a 
transiting exoplanet. Using Spitzer observations of 7 transits 
of HD 189733b, we have placed upper limits on the planet's 
oblateness. The observed absence of lightcurve anomalies 
constrains the sky-projected oblateness and obliquity at each 
epoch. The collection of light curves spanning 1 .6 yr con- 
strains the true oblateness and obliquity, because spin preces- 
sion would have produced transit depth variations in contra- 
diction of the data. In both analyses an oblateness as large as 
that of Saturn could be ruled out with >95% confidence. 

7 |http : //kepler . nasa ■ gov/] see also Borucki et al. (2009). 



The resulting upper limits on the oblateness are physically 
meaningful, in the sense that they represent a physically pos- 
sible degree of oblateness comparable to that of giant planets 
in the Solar system. However, for HD 189733b one would 
naturally expect a slower rotation rate and a smaller oblate- 
ness than those of Saturn or Jupiter, because tidal effects are 
expected to have slowed down the planetary rotation until it 
was synchronized with the orbital motion. Assuming that the 
planet is indeed spin-orbit synchronized with zero obliquity, 
we were able to place an upper bound on the quadrupole mo- 
ment of the planet's density distribution: < 6.8 x 10~ 2 with 
95% confidence. However this is a weak upper bound, in the 
sense that the theoretically expected value is at least an order 
of magnitude smaller. 

The planet HD 189733b was chosen for this analysis be- 
cause of the favorable stellar brightness and transit depth, as 
well as the large collection of high-quality data that are cur- 
rently available for this system. However, this system suf- 
fers from the drawback that the theoretically expected oblate- 
ness is smaller than detection thresholds. The transiting 
planets HD 17156b (Barbieri et al. 2007) and HD 80606b 
(Moutou et al. 2009, Fossey et al. 2009, Garcia-Melendo & 
McCullough 2009) are attractive targets because they have 
much longer periods and higher orbital eccentricities than 
HD 189733b; although tidal evolution is expected to be im- 
portant in both of those cases, it may have resulted in "pseu- 
dosynchronization" (Hut 1981) or some other state besides 
spin-orbit synchronization. 

Furthermore, the Kepler mission 7 and other surveys for 
transiting planets should soon find planets with longer orbital 
periods and larger orbital distances, for which there is no rea- 
son to expect spin-orbit synchronization. There should be a 
"sweet spot" in orbital distance, far enough that the planet 
might be expected to rotate as rapidly as Jupiter or Saturn, yet 
close enough for the stellar gravitational field to cause rela- 
tively rapid spin precession that would be manifested as tran- 
sit depth variations. A more general analysis of precession- 
induced transit depth variations seems warranted. 



We are very grateful to Heather Knutson, Eric Agol, and 
their collaborators, for organizing the Spitzer observations. 
We are indebted to Dan Fabrycky for pointing out that spin 
precession can play an important role in the analysis. We 
thank Darin Ragozzine and the referee, Jason Barnes, for pro- 
viding timely and detailed comments on the manuscript. We 
also benefited from discussions with Sara Seager and Saul 
Rappaport. 



APPENDIX 

A NUMERICAL ROUTINE FOR FAST, EFFICIENT CALCULATION 
OF TRANSIT LIGHT CURVES FOR OBLATE EXOPLANETS 

In this appendix we describe a method for fast and stable calculations of light curves of ellipsoidal exoplanets, for which 
the sky projection is an ellipse. The development of this algorithm was crucial for our analysis, to make the MCMC method 
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FIG. 6. — Constraints on the true oblateness and obliquity of HD 189733b, based on Spitzer observations of 7 transits. The solid blue curves bound the regions 
containing 68% or 95% of the MCMC samples in the oblateness-obliquity (f-9) plane. The black points are a random subsample of the full Markov chain, shown 
to illustrate the posterior probability density. The short-dashed green contou rs in dicate the amplitude of transit depth variations. The long-dashed red contours 
show the expected spin-precession period, according to the formalism of § 12.31 The hatched region was not included in our analysis because the estimated 
precession period is > 10 yr. 

computationally tractable. Proper convergence of these posterior distributions required many millions of executions of the light 
curve calculation. Seager & Hui (2002) and Barnes & Fortney (2003) both calculated light curves of oblate planets, but for those 
authors the computation time was probably of secondary importance. 

Figure|7]illustrates the basic geometry of the calculation. An elliptical shadow with major and minor axes a and b obscures a 
portion of the circular stellar disk. The projected-oblateness f± is defined as (a-b)/a. Distances are expressed in units of the 
stellar radius, R+. The semimajor axis of the ellipse is inclined by an angle a from the line connecting the centers of the ellipse 
and circle. The angle a is generally not the same as 9±, but the angles are equal for a central transit. The centers of the star and 
planet are separated by a distance x. We assume the stellar disk has a brightness profile 7*(r, 9) including, for example, a radial 
limb-darkening profile. The fractional flux deficit, F(x;a,b, a,h), due to the obscuring ellipse is 



F(x;a,b,a,L)=-^ I h(r,8) r dr d6 (Al) 
lenc 



F 



where Fq is the total unobscured flux and the integral is performed over the region bounded by the intersection of the ellipse and 
circle, denoted as £ n C. We consider a radial brightness profile given by a quadratic limb-darkening law, 



L(r,9) = I+(r,u l ,u 2 ) 



= /*(!) 



-Ml (l-Vl-r 2 ^-M 2 (l-Vl-r 2 y 



(A2) 



whereMj and«2 are the limb darkening parameters. With this law,Fo//+(l) = 7r(l-l/3«i-l/6M2)- Wedenotethe quadratic-profile 
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fractional flux deficit by the form F(x;a,b,a,u\,ii2). For a = b = R p /R±, F(x;a, u\,ui) = F(x;a,a,a,u\,U2) has a closed-form 
analytic solution in terms of elliptic integrals (Mandel & Agol 2002). The problem is more complicated for arbitrary a, b, a, u\ 
and «2 since the vertices of the intersection region £ DC cannot always be determined analytically, and the integral cannot be 
done in closed form. 

There are standard routines for computing the intersections of two ellipses that would allow us to define the intersection region 
£ DC (see, e.g., Hill 1994). We found, however, that these algorithms are not numerically stable over the entire region of parameter 
space of interest. We chose to sacrifice runtime to ensure stability by finding the intersections of a polygonal representation of the 
obscuring ellipse and the disk, as follows. First, we analytically find the coordinates in the x—y plane of points on the boundary of 
the ellipse at N uniformly selected angles (from to 27r). Next, these points in the plane are connected by line segments. Finally, 
the routine accumulates approximate intersections of the ellipse and circle by determining the intersections of each line segment 
and the circle. We use N w 200 to ensure adequate precision. 

With the points of intersection determined, we move on to the calculation of the integral in Eqn. (IA1I) . Before discussing the 
general problem , we first consider the case of uniform brightness (ui = U2 = 0) and solve for F(x;a,b,a,Q,Q). Here, solving the 
integral in Eqn. (IAU is equivalent to finding the area of the region £ DC. This area may be calculated analytically by using the 
formula for the area of an elliptical chord A{9\ , 9 2 , a, b): 

A(0i ,e 2 ,a,b) = ab[(9 1 - 8 2 ) - sin(0, - 2 )] /2 (A3) 

where the angles 9\ 2 are measured relative to the semi-major axis. To solve for F(x;a,b,a,Q,0) we add the area of the elliptical 
chord (defined by the lines connecting the intersection points 8 and the curve bounding the ellipse that is internal to the stellar 
disk) to the area of the circular chord (defined by the complement of this elliptical chord and £ flC). 
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FIG. 7. — Geometry of the transit of an ellipsoidal planet across a circular stellar disk. The deficit in flux due to the obscuring ellipse with semi-major axis a, 
semi-minor axis b and orientation a a distance x from the center of the star may be found by integrating the brightness profile of the star over the intersection 
region of the ellipse S and circle C. The ellipse and circle intersect at points x\2- A closed-form solution for the flux deficit exists for the inscribed circle C of 
radius b. Only the intersection region less this circle, £ /C C\C needs to be integrated numerically. 

Next, consider the integral given by Eqn. (lAlb for nonzero u\ and m. Numerical integration must be performed. We could 
perform two-dimensional numerical integration over £ (~l C [as was done by Seager & Hui (2002)] or numerically integrate in 
the radial direction, calculating the intersections of circle and ellipse at each integration step [as was done by Barnes & Fortney 
(2003)]. For our application we chose the former, although because of the difficulty in specifying the integration region, we found 
standard deterministic integration techniques to be too slow. 

Instead, we chose to use Monte Carlo integration [see, e.g., Press et al. (2007)]. With Monte Carlo integration, the integral 
J A f(x,y) dx dy over the region A is calculated by sampling the function f(x,y) at jV uniformly distributed random points in a 
region 1Z that covers A. When a sample (x,y) is drawn from 1Z that is not in A, we set f(x,y) = 0. According to the fundamental 
theorem of Monte Carlo integration, 



/ f( x >y) dx dy ~ Area[7\L] 
J A 

8 In general, a circle and an ellipse may intersect at up to four distinct locations; for an exoplanet with small oblateness and R p Sj up i lcr orbiting a star with 
i?* Rq , only two intersections (at most) are expected. 



</)± 



Var(/) 
N 



(A4) 
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The means are calculated over the N randomly sampled points (Press et al. 2007). In general, Monte Carlo integration is reserved 
for integrals of high dimensionality due to its slow convergence (1/ \/N) compared to more traditional methods when calculat- 
ing two-dimensional integrals. For our purposes, Monte Carlo integration reduces the difficulty associated with bounding the 
intersection region £ DC. An improvement in convergence (~ 1 /N) is obtained by sampling 1Z "quasi-randomly" as opposed to 
randomly, using a low-discrepancy random sequence called the Sobol' sequence (Sobol' & Shukhman 1995, Press et al. 2007). 
The Sobol' sequence is a sequence of uniformly distributed values on the unit interval, s, for i > 0, such that a given sequence 
member at index / is maximally distant from all previous samples i < I. 

To simplify the computational effort required for the integral of Eqn. (IAU when using the quadratic limb-darkened profile 
[Eqn. (IA2H . we rewrite the integral as 



, Ml «2 

F(x;a,b,a,ui,U2) x tt I 1 — — 

3 6 



l- Ml (l-Vl-r 2 )- M2 (l-Vl-r 2 ) 



2 

r dr dO 



£nc 

= (1 — Mi — W2) ( / r dr dO 

Unc 



+ 



r dr dO 



(mi + 2u 2 ) Vl-r 2 - u 2 (l- r 2 ) 
(I-M1-M2) x F(x;a,b, a, 0,0) 

I i ,(r;ui,U2) r dr d9 (A5) 

snc 

where F(jc;a,fe,a,0,0) is the uniform-brightness solution (calculated as described above) and U{r;ui,U2) = (mi+2m2)V1 — r 2 — 
112 ( 1 - r 2 ) describes the nonuniform component of the brightening profile. IJr = 1; U\ , U2) = at the stellar limb. The remaining 
unknown integral is over a function that has low variance over any bounding region 1Z covering the limb, making the absolute 
error in the Monte Carlo integration smaller (for a fixed AO as compared to integrating in full. 

Note that the only portion that requires numerical integration is the intersection region excluding the circle C of radius b 
inscribed in £ [see Fig. (|7)]. We denote this smaller region by £/CnC. The integral over C is analytic, as per Mandel & 
Agol (2002), and is denoted F(x;b, 111,112). Thus Eqn. ( IA5b becomes 



F{x;a,b,a,u\,U2) x tt (\ 



U\ 112 " 

(1 - Ui - u 2 ) X [F(x;a,b, a, 0, 0) - F(x; b,0,0)] 
+F(x;b,ii\,U2) 



+ / /*('"> Mi, M2) r dr d8 '. (A6) 

Jsr\C/c 

The only nontrivial component is the final term, which we denote I and for which we use quasi-Monte Carlo integration. For 
slightly oblate exoplanets or weakly limb-darke ned brightness profiles, the contribution of this integral in the total flux deficit is 
small compared to the remaining terms in Eqn. ( 1A6I >. The absolute contribution of this term to the flux deficit is bounded by 

HV-T-T) J£ i-i/3Mi-i/6» 2 x ^- fc ) 

=abf± — U ; +2U2 / . ( A7) 

1 — 1 /3m 1— 1 /6M2 

For comparison, the size of the flux deficit in the absence of limb-darkening is ab. 

To calculate X via quasi-Monte Carlo integration, we select a covering region, 7Z D £/C DC, that efficiently bounds the inte- 
gration region and that can be easily sampled with the Monte Carlo technique (£/C n C cannot). Here, easily sampled means that 
we may take a two dimensional Sobol' sequence, uniform in [0, 1] x [0, 1], and analytically transform it such that the transformed 
sequence uniformly samples our chosen covering region. A uniform sampling on the unit square may be mapped to a uniform 
sampling over a region bounded by an elliptical annular sector (see AppendixlBlfor details). An elliptical annular sector is the re- 
gion bounded between two concentric ellipses (i.e., two ellipses with coincident centers and equal axis ratios and position angles) 
with semi-major axes a\ < 02 and by the rays (emanating from the common center) at angles 9\, O2 relative to the semi-major 
axis. The formulae for the specific values of 9\ , 62, a \ and «2 as functions of the parameters x, a, b, and a for use in the integration 
routine are chosen from one of the following cases [with reference to Fig. (8j]: 

(7) One or fewer intersections between circle and ellipse. 

(a) x > 1; The obscuring ellipse is external to the circle, X = 0. 

(b) x < 1; The obscuring ellipse is properly contained in the circle; 
9\ = 0, 62 = 2tt, a\ = b, 02 = a. 

(II) Two points of intersection between circle and ellipse. 
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The circle and ellipse intersect at ^-coordinates X12 at angles 9\ 2 from the x-axis. 

(a) x > 1; For each intersection x-coordinate x' = x-„ angle 9' = 9[ 
and final solution angle 6 = 8,: 

(1) \x\ > 1 1 jx' I ; The line connecting the point of intersection and the ellipse 
center [at (x,0)] intersects the circle exactly once: 9 = 9' -a. 

(2) \x\ < |l/x'|;The line connecting the point of intersection and the ellipse center 
[at (x,0)] intersects the circle twice: use 9 = 9- a where the angle 

9 = tan -1 1 / Vl-x 2 defines the line segment that connects 

the ellipse center and the point of tangency on the circle. 
If a' is the semi-major axis of the concentric ellipse "kissing" the circle 9 then 
ci\ = max(b,a'), = a. 

(b) x < 1; Use the angle 9b defined by the intersection points of the inscribed circle 
with radius b (C) and the stellar disk (C): 9\ 2 = ±#z,-a, a\ = b,a^ = a. 

By choosing the elliptical annular sector as defined above as the bounding region for the the desired integration region {& /CC\C), 
we ensure that the integration region covers approximately 50% of the bounding region regardless of the values of x, a, b, and a. 
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FIG. 8. — Quasi-Monte Carlo integration of the nontrivial component of the total flux deficit for the stellar transit of an oblate planet. These figure s sh ow the 
transit phases and quasi-Monte Carlo integration regions described in Appendix [A] which are needed to evaluate the nontrivial integral X [see Eqn. jA6H . The 
labels in the upper left hand comers of each figure correspond to those in the text. In each figure, the inner ellipse with semi-major axis a \ gives the inner boundary 
of the elliptical annular sector. The red rays emanating from the ellipse indicate the angular extent of the sector (with angles Q\i relative to the semi-major axis). 
The blue or black points are 1000 uniformly distributed Sobol' points in the elliptical annular sector. The blue points are those which fall in the integration region 

s/cnc. 



To give some idea of the computation time, we consider the case x = 1, a = 0.155, b = 0.148, u\ = 0.2, 112 = 0.3, and a = 0.5, 
corresponding to an oblate (f± = 0.05) version of HD 189733b. Using a C++ implementation of our algorithm on a 2.6 GHz Intel 
Core 2 Duo MacBook Pro, it takes 0.5 ms to compute F with a precision of 1 ppm. For comparison, it takes 0.004 ms to execute 
the same computation for a spherical planet, using the Fortran implementation of the code by Mandel & Agol (2002). 

UNIFORM SAMPLING OF AN ELLIPTICAL ANNULAR SECTOR 
Elliptical annular sector 

A point (x,y) is inside the elliptical annular sector centered at (0,0) with semi-major axis in the x-direction, axis ratio e, inner 
radius a\, outer radius 02, and sector angles 8\2 if the all of following conditions are satisfied: 

x 2 y 2 
a\ (ea{f 

9 "Kissing" ellipses intersect at exactly one point. The "kissing" ellipse's semi-major axis a! for this application is found with a linear search. 
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2 2 

(2) _ + -2- 5 <l (B2) 

«2 (M2) 

(3) The line connecting (0,0) to (x,y) is at an angle 9 relative to the x-axis such that 62 > > 9\ (B3) 
An elliptical annular sector is illustrated in Fig. 




Uniform sampling 

Let (m,v) be a uniform sample of the unit square [0, 1] x [0, 1]. Then («', v') is a uniform sample of the elliptical annular sector 
where 

u' = a 2 rcos(8) (B4) 



v'= — sin((9) (B5) 

e 



and 



r=\J(\ -u)a\ + u (B6) 
6» = ( 1 - v) tan" 1 (e tan 6 1 ) + v tan" 1 (e tan 9 2 ) . (B7) 
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